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"<n : 

The Monte Carlo (MC) estimates of thermal averages are usually functions of system control 
£SJ ' parameters A, such as temperature, volume, interaction couplings, etc. Given the MC average at 

a set of prescribed control parameters Ao, the problem of analytic continuation of the MC data to 

3 ; 

h—j ' A-values in the neighborhood of Ao is considered in both classic and quantum domains. The key 

result is the theorem that links the differential properties of thermal averages to the higher-order 
i 1 1 cumulants. The theorem and analytic continuation formulas expressed via higher-order cumulants 

O |. are numerically tested on the classical Lennard- Jones cluster system of N = 13, 55, and 147 neon 

g , particles. 

^ ■ 

I I. INTRODUCTION 

o ■ 

c/2 , To obtain Monte Carlo (MC) estimates of thermal averages, one often has to run MC codes multiple times in order 
to get the corresponding results at different values of system parameters, such as temperature, volume, magnetic or 
i— r electric fields, etc., or at different values of intcrparticle interaction constants. To avoid these extra time-consuming 
y—{ ' runs, several highly effective strategics for sampling phase space of the system and/or extracting maximum relevant 
, information from the MC data acquired have been suggested. Among them, these are histogram methods by Ferrenberg 
Q , the multicanonical method by Berg and Neuhaus J] , 



and Swendsen 



the Wang-Landau method [5| and others. 



Thus, in the first method we record a histogram n(E) of how many times each particular value of the energy E 
is generated in MC simulation at a particular temperature To. Then, using this energy distribution histogram one 
can in principle recalculate the corresponding energy distribution and thermal averages at an arbitrary temperature 
J> ' T. The multicanonical method is based on the idea of using inverse of the density of states, l/p(E), instead of 

•1— I ' 

the Boltzmann factor exp(— j3E), where /3 = l/(fcgT) and fee is the Boltzmann constant, for sampling the energy 

5_i , 

C$ ■ states in the metastable-unstable region of the canonical ensemble. This results in a histogram in which all energies 
are sampled equally. The obvious problem with the direct implementation of this idea is that we do not know the 
density of states p(E). However, using a sequence of approximations, Berg and Neuhaus were able to demonstrate 
that p{E) can be found, even in a particularly impressive case of a first-order phase transition [4(. On the other hand, 
the Wang-Landau algorithm allows to estimate the density of states p(E) directly instead of trying to estimate it 
from the probability distribution obtained at To and then with density of states one can easily calculate the partition 
function, free energy, etc. at an arbitrary temperature T. Moreover, the original Wang-Landau algorithm |5| proposed 
for MC sampling in spin lattice systems has been generalized to the off-lattice systems, such as continuum (fluid) 



models [6), [7[ , polymer films 



By a suitable reformulation of the problem the Wang-Landau sampling scheme can 



be advantageously employed in quantum systems as well 

However, it should be noticed that the applicability of these methods is severely restricted by the presence of 
statistical errors in the data, n(E) or p{E), generated in MC sampling. If errors in the data become comparable or 
bigger than the true values in the energy range near the average potential energy at temperature T, then the above 
methods fail to continue the corresponding thermal averages from To to T temperatures. 

In statistical physics, quantum mechanics, quantum held theory, and many other fields of modern theoretical physics, 
where some kind of averaging of a generalized exponential function is present, we constantly face the so-called cumulant 
expansions 



1CH12|. Although cumulants (semi-invariants) have been known long before in mathematical statistics 
and probability theory, it was Kubo who first convincingly demonstrated in his pioneering work [13I ] how the concept 
of cumulants can be widely applied to various problems of quantum mechanics and statistical physics. For example, it 
has been successfully applied to Ursell-Mayer expansion for classical and quantum gases [14j that is usually obtained 
by much longer diagram considerations, to perturbation series in quantum mechanics and random perturbations 
in dynamical systems, to relaxation functions in irreversible processes [13I [l^j], etc. In spite of such a diversity of 
applications discovered so far, cumulants and their properties have not been, to the best of our knowledge, exploited 
before in the context of analytic continuation problems both in classical and quantum domains. In his original 
paper, Kubo greatly generalized the classical concept of cumulant expansion to the case of non-commuting operator 

I L 

generating algebra and applied relations found in this algebra to a diverse set of physical problems. In 16], a number 
of useful algebraic and geometric properties of cumulant expansions have been summarized and applied to generate 
cumulant Faddeev-like equations and to establish a method of increments for excited states. Very recently, cumulant 

However, application of 



expansion techniques have been successfully applied to the Fourier path integrals 
cumulants to the problem of analytic continuation requires the knowledge of differential properties of cumulants. 
These properties have not been established before. Therefore, one of the goals of the present paper is to derive those 
properties of cumulants that are of paramount importance for analytic continuation applications. 

In this work, we present a new, more robust cumulant expansion method which allows to analytically continue 
thermal averages in the neighborhood of To. In particular, we prove the key Theorem II V. 1 1 which relates the derivative 
of the fcth order cumulant with respect to the inverse temperature /3 to a (k + l)th order cumulant. Based on this 
theorem, one can develop an asymptotic expansion in the neighborhood of To m terms of higher order cumulants. 
Also, some numerical results for the heat capacity of the Lennard- Jones (LJ) clusters illustrating the cumulant's 
derivative theorem and the quality of the derived expansions are presented. One of the advantages of the proposed 
method is that the analytic continuation can be implemented directly without need of separate calculating n(E) or 
p(E) distributions. Moreover, given the cumulants, if necessary one can easily express these distributions in terms of 
cumulants. 

The rest of the paper is organized as follows. The energy representation for the phase-space probability density func- 
tion (pdf) is introduced in Section UH In the energy representation, all degrees of freedom irrelevant to thermodynamic 
equilibrium are integrated out. Next Section considers how partition function, entropy, and statistical temperature 
can be expressed via kinetic and potential energy pdfs. The thermodynamic energy, heat capacity and its differential 



3 

properties are given in terms of cumulant expansions in Section llVl Here, the key Theorem IIV.1I about cumulant's 
derivative is formulated. In Section |Vj the cumulant expansions and the corresponding analytic continuation formula 
for the energy pdf are analyzed. The statistical or microcanonical temperature in terms of cumulants is analyzed in 
Section IVII Further generalizations to the multi-parameter classic and quantum systems are developed in Sections 
IVIII and lVIIII Here, Theorem lVII.il a multi-parameter generalization of the Theorem HV.ll is formulated. Numerical 
results are discussed in Section [IXI Concluding remarks are in Section [X] Finally, technical details about cumulants, 
proofs of Theorems IIV.ll and IVII. II can be found in Appendices [AllCl 



II. PHASE SPACE TO ENERGY SPACE MAPPING 



Let us consider a classical system, the phase space of which is described by a set f2 of 3 AT canonically conjugate 
coordinate r = (ri, . . . , rjv) and 37V momentum p = (pi, . . . , pjv) variables. On the phase space we define the pdf 
p(Vt) that fully describes a thermal equilibrium state such that the thermal average of an observable 0(fi) can be 
calculated as 

(o(n)) = J d m no{rt) P {n) (i) 

where d 6N il = d 3N rd 3N p/w]y and is an appropriate weight factor, sometimes called Gibbs factor, which makes 
the classical averaging as close as possible to the quantum-mechanical one. Further, we assume that the pdf depends 
on K functions H(fl) = (Hi(£l), . . . , Hk(Q)) and K control parameters A = (Ai, . . . , A^) so that its Sl-dependence 
can be represented as p(Q) = p(H($V),\) For the observable we assume a similar structure, 0(£1) = 0(H(n)), to be 
valid. If we define the ^-dimensional energy pdf as 

K 

P (E, A) = / d 6N n Yl S(E k - H k {Q))p{H(Q), A) (2) 

fc=i 

where E = (E\, . . . , Ek ), then, the thermal averages can be be evaluated by integration in the energy space only 

(O(E)) = J d K EO(E)p(E,X) (3) 

The original problem of integration in a 6./V-dimensional phase space is, thus, reduced to a i^-dimensional integration 
in the energy space. In Eq. (|2]), we effectively integrated out all extra degrees of freedom not important in the 
equilibrium state. 

The problem that now can be formulated is how to calculate this energy pdf at a fixed set of parameters. The 
energy pdf at an arbitrary set of parameters A can be calculated, at least in principle, if it is known at a fixed one, 
Ao, as follows. Usually, the pdf in the phase space is known up to the normalization factor or partition function. We 
assume a generic form for 

p = exp(-\- H(Cl)) /Z(\) (4) 
where • denotes the scalar product and the partition function 

Z(A) = J d 6N n cxp (-A • H(Q)) (5) 



4 



It assumes the existence of the partition function. 



20 



2l\ . for which there is no need a priori to know 



Using Metropolis et al. importance sampling algorithm 
the partition function, one can evaluate the integral over Q by the Markov chain MC simulation method and get 
a statistical estimate for the energy pdf p{E 1 Ao) at Ao parameters. To this end, one can use a proper asymptotic 
representation for the ^-functions in the integrand of Eq. ([2]) 

6{E k -H k ) = - hm . (6) 

7T r s-oo hj k - ti k 

Having obtained the energy pdf at a fixed value Ao, one can easily recalculate the pdf at an arbitrary value A using 
the formula 

p(E, A) = «»P(-(A-y-g)p(S,Ao) (?) 
Z(X, A ) 



where the normalization factor 



Z(A,A ) = f d K Eexp(-(\-\ )-E)p(E,\ ) 
Z{\) 



(8) 



Z(\ ) 

is the ratio of partition functions calculated at A and Ao parameters. Thus, the thermal averages at arbitrary values 
of A can be evaluated with the help of Eqs. © and (J7J. Examples of such calculations in the system of particles 
interacting via L J potential will be presented in Section IIXI 



III. CLASSICAL PARTITION FUNCTION, DENSITY OF STATES AND STATISTICAL 

TEMPERATURE 

At first sight, the best one can get from Eq. is the ratio of partition functions, not an absolute value at a particular 
set of parameters. However, at A = the partition function takes an especially simple value: Z(0) = J d 6JV Sl = Vq/wn, 
where Vq is the total available volume of the phase space. In classical, non-relativistic physics Vn is, in principle, infinite 
due to possible infinite particle's momentum values. However, it is well known [22| that kinetic energy contribution 
to the partition function, Z/c(f3), where /3 is the inverse temperature, can be evaluated exactly and calculation of the 
partition function, Zx,+v{\) = Z)c{fi)Zy(X), is, thus, reduced to the computation of the configuration integral, iJy(A), 
defined as in Eq. ([5]), where Vt — > r is a position in the configuration space and H(£l) is replaced by the potential 
energy V(r). 

Therefore, the configuration integral 

Zv(P) = V N , (9) 

where V is the spatial volume of the system. From Eq. ([8]), we find that 

7 ,^_^ N Jd K Ee X p(-(\-\ )-E)p v (E,\ ) 

ZvW ~ V j d K E exp(A .£)M£,Ao) (10) 



Here, pv(E, Ao) is the potential energy pdf. Notice that at A = Ao, the integral in the numerator is equal to one due 
to normalization of the pdf, whereas at A = 0, Eq. (fpQ|) is reduced to ([9]). 

As an example, let us consider the case of a single parameter Aoi = /3o = 1/(&b To), where To is an absolute 
temperature, Vi(r) = V(r) is the potential energy. We wish to calculate the density of states 



Pic+v 



(E) = J d 6N nS(E-IC{p)-V(r)) (11) 



In contrast to the partition function, the calculation of the density of states cannot be factorized in separate integrals 
over momenta and coordinates. To overcome this difficulty, it is useful first to map the phase space to the two- 
dimensional energy space as suggested by Eq. @, namely, to define the 2D pdf in a factorized form 

Picy(Ei,E 2 ,f3 )= p^{E 1 ,p o )p v {E 2 ,0o), (12) 

exp(-A)/C(p)) 



p K (E 1 ,/3 )= /d 3JV p5(£i-£(p))- 



Z K (Po) 



M E„M= /^r^-V(r))5Eb^p 



23], Ch. 3. 



Similar to the kinetic energy partition function, the pdf px; can be calculated exactly; see, e.g., 

With the 2D energy pdf (fl2]) . Eq. (fTT|) can be rewritten as a convolution of the kinetic and potential energy pdfs 

p K+v {E) = Z K+v (/3 )exp(/3 S) 

x J dE 2 p K (E-E 2 ,f3 )p v (E 2 ,f3 ) (13) 

Substituting explicit expressions for Z/c+v ancl Pk into f) 1 3 j) . one obtains 

V N (2nm)*? 
PK+V{E) = ~^7F3N 



2 
E 

3N 



r(^) 

J dE 2 (E-E 2 )^- 1 eMPoE 2 )pv(E 2 ,(3o) 

■OO 

OO 

J dE 2 exp(p E 2 )p v (E 2 ,(3 ) 



(14) 



24] , expressed in terms of the potential energy pdf. In 



where m is particle's mass and T(^j-) the Gamma function 
the limiting case of zero interaction, V(r) = 0, pv-fo(E) —> 5(E) is reduced to a (S-function and the density of states 
Eq. dHD is defined only by the kinetic energy contribution, p K {E) = 6(E)V N (2nm) 3N / 2 E 3N / 2 - 1 /T(3N/2), where 
0(E) is the Heaviside step function. It is easy to check this result directly from Eq. (fTT|) . 

Moreover, from Eq. (fT4|) for the density of states we find the corresponding expressions for the entropy, S/c+v(E) = 
ks In pk.+v (E) , and the statistical temperature 

fdS^MEiy 1 9l (E) 

Tk+v{e) = {-^e—J = MFW' (15) 

E 

9 P (E) = J dE 2 (E-E 2 )^-Pexp((3 E 2 )p v (E 2 ,(3o) 

— CO 

OO 

/3N _ 
dE x E^ p exp(-f3 E 1 )p v (E~E 1 ,(3 ) (16) 



in terms of the potential energy pdf. If V = 0, then from Eq. (fT5]) one gets for the kinetic energy temperature 



6{E)E 



Tk{E) = - (17) 



Here, although formally pic = and Tjc does not exist at E < 0, we put T/c — by a continuity at negative energies. 
Then, we can measure the potential energy contribution to the temperature as the difference Ty = Tjc+v — Tjc- 



IV. THERMODYNAMIC ENERGY, HEAT CAPACITY AND CUMULANT EXPANSION 

The thermodynamic average of the energy written via the first moment or the first cumulant term reads as 

Utf) = |^+Mci(/3), (18) 
/x cl (/9) = M i(/3) = (E) ME ,p) = ( V (r)) Pv (r, P) (19) 

where the first and second terms are, respectively, due to the kinetic and potential energy contributions. Here, (. . .) 
denote averages over the potential energy pdfs either in the energy or coordinate spaces. 

Differentiating the fcth-order moment /ifc(/3) = , B ^ with respect to (3, one obtains from the definition 

— = -fik+i + MfeMi (20) 
The following theorem establishes a similar relationship between the derivative of the fcth and (k + l)th cumulants. 

Theorem IV. 1 (univariate): Let fj, c k(/3) be the kth-order classical cumulant. Then, 

d^ck ,„-.. 
-7]r = -Mc(fc+i)- (21) 



The general definition of the fcth-order cumulant is given elsewhere, see, e.g., [13|, and Appendix lAl For proof of 
Theorem IIV. II we refer to Appendix [Bl 

Using the relationship between the derivatives d/dT = —ksfid/ 'd/3 and Eq. pTj) . one can easily find derivatives of 
energy (|18p with respect to temperature expressed in terms of the second and higher-order cumulants. For example, 
for the first two derivatives we have 

= d ~l^ = fc|/?3 h2/ic2(/3) + (22) 

From (|22p one immediately derives that the extrema points of the heat capacity curve are defined by equation 

2 Mc2 03) = jfyios 09) (23) 
Using Taylor's series expansion near a fixed value /3q and the equation for derivatives 



d k ' 



df3 k 



(-l)>c( P +fc), P,k= 1,2,3,... (24) 



that follows from Theorem lIV.il Eq. (|23l) can be rewritten as 



E [PoVc(k+3)(M ~ (k + lK( fe+ 2)(A))] ( ^ = Mc 2 (/3o) (25) 

k=0 

where Aft = ft — /3q. Solving this algebraic equation truncated at some maximum power k — k max for Aft, one obtains 
a converged root that can locate an extremum position. 
Similarly, for the heat capacity one obtains the expansion 



^ = |iV + ft 2 [ Mc2 (A)) - ^ c3 (/3o)A/3 
kb 2 

+ ^ c4 (/3 )A/3 2 - i Mc5 (/3 )A/3 3 4 
2 o 

Some numerical applications of these equations will be given in Section 



(26) 



V. THE ENERGY PDF AND CUMULANT EXPANSION 

Using the Fourier integral representation for ^-function, we obtain an integral representation for 

oo 

p v (E,T) = ±- J drexp(irE)(exp(-iTV(r))) pv{rT) (27) 

— OC 

Formally, averaging the exponential function can be rewritten in terms of cumulant's expansion as 

(exp(-irV(r))) pv(r T) = exp (S c (r, ft)) , (28) 

k=l 

where fi c k is the fcth cumulant. On the other hand, with the help of Taylor's series expansion and Theorem lIV.ll one 
obtains the following cumulant expansion for 

„ , fl . lVl f^(ir) fc dV cl (/3) 

» cl (ft + lT) = ^ — -g-, 

= E-^-^m (29) 

k=0 

Therefore, from (12511 and (l2l?l) one derives alternative integral representations for 

E c {r,ft) = -i [ dr'^ift + iT 1 ) (30) 

and 



(i 



oo 

Pv(E,T) = -±- J dTexp(iTE-i dT'fidift + ir')) 



(31) 



Let us consider truncated cumulant expansions defined as 



)(r,ft)= f^L^tf, ck (ft) (32) 

k=l 



1 we easily get a <5-function like pdf 
j#> =S(E- Mcl 09)), (33) 



while at fc mQX = 2 integration over r yields the Gaussian function 

(£- Mcl (/3)) 5 



(2) 1 

Pv = — z exp 



(34) 



2^ c2 (/3) 

The next case of k max — 3 is more complicated; the pdf is reduced to an Airy function. First, we have to shift the 
integration variable to the complex plane r 1 = r + i^ C 2/ 1 Hc3 in order to cancel the quadratic term in the exponent 



Py } = J drcxp |ir(E - (i cl ) - y^ c2 + iyAi c3 | 



00 + i/i c 2/pc3 

2 1/3 exp(0) f , f 

— — / exp ^ 



-00+1^2/^3 

where 



27r|Mc3|1/3 y dr'expjW + i-^ (35) 



(£/ - ^cl) + ^-g-. 



m c2 



1/3 
Mc3 



(36) 



Let us consider a closed rectangular path C = U»=i Q m the com plex plane r [see Fig. [TJ consisting of two 
finite horizontal G\ and C3, and two vertical C2 and C4 segments. The horizontal segments are defined as 
C\ — {t : — t < Rer < t, Imr = ^2/^03} and C3 = {r : — t < Rer < t, Imr = 0}, where < t < 00 is a pa- 
rameter fixing the ends of segments, while the vertical segments C24 = {t : Rer = ±t, < Imr < ^2/^3} connect 
the corresponding ends of the horizontal ones. Here, we have assumed that C± lies in the upper plane, i.e., [i C 3 > 0; 
the negative case can be considered similarly. According to the Cauchy theorem, contour integral of a holomorphic 
function along a closed path is zero. Thus, we have § c dr ■■■ = Yli=i Ic dr ■ ■ ■ = 0. It is easy to check that in the 
limit t — > 00 contributions from the vertical segments go to zero and, therefore, integral along the real axis and the 
one taken along the path shifted into the complex plane turn out to be the same 

00 + iA'c2/A'c3 



dr exp |^icjr + i— j = J dr exp ^kjt + i- 

and (|35p in terms of Airy function takes the form 

( 3) 2 1 / 3 exp(</>) 

p v = ~\ iI7T~ Al ( w ) 37 

lMc3| 1/J 

(2) . . ..... (3) 

Contrary to the p v pdf, which is a symmetric distribution with respect to [i c i, exhibits an asymmetric behavior. 
The Airy function shows qualitatively different behaviors: an oscillatory one at u) < 0, while an exponential decay at 

(3) 

lu > 0. However, there is no guarantee that always p v > because of possible oscillations at oj < 0. 



One can see that the pdf py 1 "^ is not reducible to elementary functions at fc max = 3. Moreover, it is expected that 
in general j5y m cannot be expressed via elementary functions at fc max > 3 as well. At fc max > 3, the asymptotic 
saddle-point (SP) approximation 
approximation, an integral 



25] can be applied in order to derive elementary working formulas. Thus, in the SP 



4 max =Re J dra(T)exp($ fcmax (r)) 



is estimated as 



Re a(i 



2n 



exp($ fcmax (r c )) 



(38) 



(39) 



where $fc max = 'wE + T EE,{ kma:c> (t, j3) is a complex phase function truncated at fc max power and the critical points r c are 
solutions of the equation $' fc (r c ) = 0. 

At fc max = 2, the SP approximation reproduces the exact result ([34]) . Let us apply Eq. (l39l) to the case of fc max = 3. 
We have the two critical points 



1T C ± 

D 



Mc3 



(40) 



2(E-fi cl ) 



Mc3 



The phase function $3 truncated at third power and its double derivative are reduced at these points to 

' . (4i) 



v ; 2 6^3 3 



$ 3 '(r c± ) = ± Mc3 V^ 
Substituting these equations into ([3T)]) . one gets 

(3SP) 



Pv 



= Rc 



1 



exp — — D tt- 



LCxp(^/ 2 ) +cxp (_^3/ 2 )]} 



3- / 3- yjj <*> 

The SP approximation (l42t is seen to have a l/4th power singularity at D = or at E S i ng = [i c \ — /i^ 2 /(2/i C 3). At 
the singularity point, where the two critical points coalesce, the standard SP approximation is not valid. 

If the pdf pv(Po) is known at a fixed value /?o, then its analytic continuation to /3 = /?o + A/3 point is given by 

exp(-Al3E) Pv (E,M 



where the normalization factor 

S(A/3,/3 ) = j dE cxp(-A/3£)p v (i;,/3 ) 
Let us calculate the normalization for the first three truncated pdfs, p^ ' . One gets consecutively 

S« = exp(-(A J 8) Atcl (A))), 



(43) 



(44) 



S (2) = exp(-(A/3K 1 (/3 ) + Alc2(/3o) - 



S {3) = exp(-(A/3)/i cl (/3 ) + 



2 

Mc2(/3 ) 



(A/3) 2 - ^M(A/3) 3 



(45) 
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where in the last line the Airy averaging has been carried out with the help of an integral formula [2n] 

/oo 
dt exp(pt)Ai (t) = exp(> 3 /3), Rep > 0. (46) 
-oo 

If the average potential energy /i c i is negative, then the normalization factor is exponentially small or large, depending 
on the sign of A/3. At A/3 > 0, S is large or it is small otherwise. 
Eqs. (|45|) suggest an anzatz for 

S(A/3,/3 ) = exp( K (A/3,/3 )), 

OO 

«(A/3,/3 ) = 5> fc (/3 )(A/3) fe (47) 
fc=i 

where the first three coefficients are defined by Eq. (|45|) as Ck(0a) = {— I) +1 fi c k(0o)/k\, k = 1,2,3. It can be easily 
checked that the same formula holds valid for the higher-order coefficients at k > 3. It follows directly from equation 

(/KRR , rf 2 ff(A/3,/3 )MA/3) 2 f dS(A(3,(3 )/d(A(3) \ 2 

rf 2 ^(A/3,/3 ) 

= d(K0y ( } 

as a result of equating the like powers in Taylor's expansions for /i c2 (A/3, /3 ) [see Eq. (|26l) ] and the second derivative 
of k(A/3,/3q) in the right-hand side. 

VI. THE STATISTICAL TEMPERATURE AND CUMULANT EXPANSION 

Making use of the integral representation Eq. ([3"Tjl . the statistical temperature in terms of cumulants can be 
rewritten as 

Tk +v (e) = T^nk, (49) 



f p (E) = I dra p (r)e X p($(r)) (50) 

where 



Op(r) = (/3o + ir)- Q -, a p = 3N/2-p + l, 

$(r) = vrE-if dT^ cl (f3 + iT') 
Jo 



fc=i 



The second line in $(r) is due to Theorem lIV.ll f p (E) can be written in an explicitly real form 

/>oo 

f p {E) =2 dr 
Jo 



eXp(g(T)) /2 cos(^(r)) (52) 

52 i ^2\ q p/ 2 



(/? 2 +T 2 r 



where 



I \ \ " Mc(2s)(A) / .... , 

7(t) 2^ ~ ' 



-r 



- (2*)! 

'(2,+ I)! 



/ \ T71 A*c(2s+l)(/3o) , , s 2S + 1 , T / K o\ 

p(r) = Et~}_^ rVTZTV ^ ^ T -a P arctan— (53) 



s=0 
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While the 'real' form (l52l might be convenient for numerical estimates of the integrals, the 'complex' representation 
(|50|) is a good starting point for developing the SP approximations. The critical points r c are defined by 

fc„ 



(-ir c ) k = E 



fe=0 



fc! 



(54) 



Solving Eq. (I54|) for t c , one gets the critical points. In the simplest approximation, truncating equation at k n 
we find 

Mc2(/3o) 

Then, with the help of (|39[) one obtains an asymptotic estimate for the temperature 

To 



(55) 



T K +v{E) 



1 



To 



k B T {E - Mci(/8o)) 



1 + 



Mc2(/?o) 

k B T {E - n cl (Po)) 



(56) 



Mc2(/3o) 

where the second line is valid if the ratio in the denominator taken by modulus is much less than one. In the 
neighborhood of /U c i, T^+y is seen to grow linearly as E increases. Notice that a similar piecewise linear interpolation 
scheme for the statistical temperature has been suggested in the statistical-temperature MC (STMC) and molecular 
dynamics (STMD) algorithms Q- 

At k max = 2, solving the quadratic equation one obtains the two critical points (|40p . Substituting the phase 
function and its double derivative (|5"Tj) taken at these points into (131?)) . where a(r) = (/3q + ir)~ ap , one gets 



f p (E) ~ReJ2 



-1 



Mc3 



exp 



V 3 



(57) 



Let us assume that E is in the neighborhood of [i c \ such that D > 0. Then, depending on the sign of /i C 3, we find 
that one term in the sum (|57|) is real and the other purely imaginary. Thus, if fi C 3 > 0, the "—"-sign term is real, 
while the "+"-sign term is imaginary. Taking into account only the real contributions to f p , one obtains 



T K+ v(E) 



Po + — - sign (Hczsj^/D 

Mc3 



(58) 



where sign (a;) = ±1 if x is positive or negative. Moreover, if 



2(£-Mci) 



and 



D 



Mc3 



Mc2 



Mc2 
Mc3 



Mc3 



sign (/x c3 ) 



E - He 

Mc2 



then (JSSJ reduces to ([561 . 
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VII. MULTIVARIATE CUMULANT EXPANSIONS 

It is straightforward to generalize the above results to the multivariate case. Thus, using Fourier integral represen- 
tation for (5-functions in the energy pdf ([2]), we obtain 



^ = (2.) 



p(o,a) 



oo / K 

= ^J---jdr 1 ...dr K e^[^E k r k+ ^ ffi^ AW,* ] (59) 

where summation runs over non-negative integers (si, . . . , s_ff) ^ (0, . . . , 0). The pdf in the phase space, p(fi, A), 
is defined by Eq. Multivariate moments can be defined as averages either over the energy or the phase space 

pdfs: fX Sl ...s K = (-^i 1 ' ' ' E s I i < ) p ( E _o = {H^ 1 ■ ■ ■ H^) p , n ^. Relationships between multivariate moments and cumu- 
lants (J-csi...sk can be established using the moment-generating function [see Appendix [U] . Truncating the cumulant 
expansion in the exponent by quadratic terms, s\ + . . . + sk < 2, the A"-dimensional Gaussian integral evaluates to 

K 

2 



P GK (E, A) = — exp 

[(27r)^deta] 1/2 



1 K 

- {E k -E k ){a- 1 ) kk ,{E k ,-E k ,) 



fc,fc'=i 



(60) 



where 



E k (\) — /io...oio...o = A 1 co...oio...o — (-^fe) p (n,A) ' 
o-fcfc'(A) = M c o...oio...oio...o = (fl^Uj) - (Hk) p rn,\) ( H k') p (n,\) ( 61 ) 

The multivariate analogue of Theorem IIV. 1 1 on cumulant' s derivative is 
Theorem VII. 1 (multivariate): Let /x CSl ... Sr ... S/< - (Ai, . . . , A r , . . . , Xk) be a K-variate cumulant. Then, 

C gx Sr " (■ • • Ar • • •) = -Mc...(» r +i)...(- ■ • A,- • ■ •)■ (62) 
Proof of this equation is similar to that done in the univariate case [see Appendix [C] . 

As an example, let us consider an application of this equation to the LJ cluster system. Let N particles interacting 
via LJ potential Vlj be put in a cubic thermostatic container of size L and volume V = £ 3 ; the system is kept at 
temperature T. The free energy of system 

F = F ldeal - fc B T In J ... J d 3N r exp (-/3V LJ (r)) (63) 

where Fideai is the free-energy of an ideal, non- interacting system and d 3N r = d 3 ri . . . d 3 r^. It is convenient to rescale 
the LJ potential to the size of container r > u = r/L separately for the repulsive and attractive parts: 



f3V LJ {r) = + PV 2 {r) 

N , N 12 N 

= **u E fen E ' 



, ,Ti — r.,- / ' — ' \ r,- — r, 

i<j=l x 1 1 J,/ i<j=l 

AiVi(u)+A 2 V2(u) (64) 
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where elj and olj are the standard LJ energy and length parameters and 

N 1 fV \ 4 

%(«) = - E^. A 2 = 4^(^j . (65) 

Here, Vlj — <?\j denotes an effective LJ volume. Notice that (i) in the rescaled variables integration over u in the 
configuration integral is carried out in a 3./V-dimcnsional hypercube of unit size and (ii) all the dependences on system 
parameters T and V are included in the dimensionless parameters A = (Ai, A2). Taking derivative of the free energy 
with respect to V at fixed T and TV, we obtain pressure 



dF\ Nk B T 



P = 

expressed in terms of the first order cumulants 



dV) T . N V 



4Al m^ 2Aa r\\ 
-^-Mcio(A) + — Mcoi(A) 



(66) 



/—/ d m uV 1 {u) exp(-AiVi(u)-A 2 F2(u)) 

Mcio(A) = (Vl) p(uX) = j , 

d 3Ar u exp(-AiVi(u) - X 2 V 2 (u)) 



f ••• J d m uV 2 (u) exp(-A 1 y 1 (u) - X 2 V 2 {u)) 

McOi(A) = (V 2 ) p(u . x) = (67) 

/•■•/ d 3Ar u exp(-AiVi(u) - \ 2 V 2 (u)) 


Note that fj, c io and /z c oi are universal functions, in the sense that their functional dependence is universal, not 
depending on the specific potential interaction parameters ulj and £lj, as well as the system parameters T and 
V . The first term in (|66p is the pressure of an ideal system with no interaction between particles. The second 
(positive) and the third (negative) terms are predominantly contributions due to repulsive and attractive interactions, 
respectively. Strictly speaking, the second (third) term contains contributions from both repulsive and attractive 
interactions, but effects due to repulsive (attractive) interactions are expected to be dominant. 

Making use of Theorem lVII.il one can expand cumulants in Taylor's series around a fixed value Ao = (A10, A20): 

Mcio(A) = Mcio(Ao) - ^ C 2o(A )AAi - // c ii(Ao)AA 2 

(Ao)(AA 1 ) 2 + Alc2 i(Ao)AA 1 AA 2 + (A )(AA 2 ) 2 + ..., 



McOi(A) = Mcoi(Ao) - Mcii(Ao)AAi - /i c02 (A )AA 



2 



+ i/i c2 i(Ao)(AAi) 2 +Ai c i 2 (A )AA 1 AA 2 + i A i c0 3(Ao)(AA 2 ) 2 + ... (68) 

where AAi = Ai — A10, AA 2 = A 2 — A 2 o and the higher-order cumulants are explicitly defined by 

A*c2o = (n 2 )-(^i) 2 , Moll = (ViV 2 ) - (Kl) (V 2 ) , 

Mcoa = (V 2 2 )-(V 2 f, h*q = (V?)-3(V?)(Vi) + 2(V 1 ) 3 , 

Mc2i = (V?V 2 )-2(V 1 V 2 )(V 1 )-(V 2 )(V 2 )+2(V 1 ) 2 (V 2 ), 

Haz = (V 1 V 2 )-2(V 1 V 2 )(V 2 ) -(V 1 ){V 2 ) + 2(V 1 )(V 2 ) 2 , 

Mcoa = <^ 2 3 )-3<y 2 2 )(y 2 )+2(l/ 2 ) 3 ,-.. (69) 
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For brevity, we dropped pdf's indication on the averaging operation. With the help of Eqs. (l66|) - (|69|) . one can 
analytically continue the state equation in the neighborhood of Ao- 
The critical point (T c , V c ) is defined by equations 

dp \ / d 2 p\ 

dVj T>N ' \dV*J 

Written in terms of cumulants, these equations for the LJ-system are equivalent to 



16A 1 /i c2 o + 16AiA 2i UcU + 4A 2 /i c o2 = N + 20\ 1 fi clo + 6A 2 ^ c0 i, 
64A?/i c30 + 96A?A 2 /i c2 i + 48AiA§Mci2 + 8Xlfi c03 = 208A^ c20 + 28A^ c0 2 + 64AiA 2 Mcii 

- 80Ai^ c i - 12A 2 ^ c0 i (71) 

Again, expanding cumulants in Taylor's series around a Ao-value, which is supposed to be close to the critical value 
A c = (Aic, A 2c ), one gets a system of two-variate polynomial equations. Once a critical solution A c of these equations 
is found, the critical volume and temperature are given by 



A 2c _ 4:£lj Aic 




V C = V LJ J-^, T c =-^-f. (72) 

V Ale k B Af c 

VIII. QUANTUM GENERALIZATIONS 

In quantum mechanics, the classical pdfs are substituted by the statistical density operators p = exp(— fiH), where 
H is the system Hamiltonian operator, so that the quantum thermodynamic average is defined by 

(Op) =tr (6 p)/tvp (73) 

where O is an observable operator. Using the path-integral representation for the density operator, a quantum system 
can be effectively mapped to a corresponding classical [polymer-type] statistical system. Based on this mapping, we 
can develop similar analytic continuation methods in quantum domain as well. Thus, the usual Feynman path integral 
expression [30( for the density matrix reads as an integral over all curves connecting the two configurations x and x': 



p(x,x';/3) = (x'| exp(-M) |x) 

/••■/ P[x(r)]exp|-^[x(r);/3]|, (74) 

x(t): x(0)=x, x(/3/i)=x' 

S[x(t);0\ = J dTH[ X (r)}= J dr | imx(r) 2 + y [x(r)] | (75) 
o o 

The symbol T> [x(t)] indicates that the integration is performed over a set of all continuous, non-diffcrcntiable [zigzag- 
type] curves x(r) : [0, (3h] — !• R d , with x(0) = x and x(/?fi) = x'; h is Planck's constant. The integer d reflects 
the dimensionality, with d — 3N for a system of N particles having a mass m and interacting via a potential V in 
3-dimensional space. 



15 



Calculating the path integral is a challenging task, which in general cannot be performed analytically. It is only 
for simple model problems, such as quadratic potentials, that an exact solution can be obtained. For more complex 
systems, the path integral has traditionally been estimated using the discretized time-slicing approximation [3l| or 



' Fourier discretization" 



32, 



33] . By introducing a change of variables to simplify the boundary conditions and 



temperature dependence: x(r) = x + (x' — x)r/ (f3h) + y(r/(/3fi,)), the reduced paths given by y, will satisfy Dirichlet 
boundary conditions, y(0) = y(l) = 0, independent of x, x', and f3. In the Fourier path representation, Cartesian 
components g/j, i = 1, . . . , d of y are expanded in a complete set of sinusoidal basis functions 

E. . . , . / -sin(/c7rw) 

a ki A k (u), A fc (u) = V? 

k=l 



kn 



(76) 



where coefficients of the Fourier expansion, {a k {\, are new functional integral variables. Or, in vector notations, we 
write y(u) — J2T=i a k^k(u), where = (a k ±, . . . , a&d). In these variables, Eq. ([74)) can be rewritten in the form 



p(x,x';/3) = lim pW(x,x';A), 

K—^oo 



P W (x,x';A) 



(K+l)d/2 



d&i 



■ dsL K cxp (-AiSf° - A 2 sf° 



K 



(x'-x) 2 +]Ta 2 



fc=i 



5'. 



(K) 



duV 



K 



x + (x - x)it + ^ a k A k (u) 



k=l 



(77) 
(78) 
(79) 



1 /2 

where Ai = 1/<t 2 , u — ((3h 2 /rn) , and A2 = (3- The <r parameter differs from the usual thermal de Broglie wavelength 
at the corresponding temperature by a factor of y/2. s[ K ^ and S% ^ are contributions to the action from the kinetic 
and potential energy operators, respectively. Observe that both s[ K ^ and S^^ do not depend on j3; all the dependence 
on j3 is separated out in the Ai^ parameters: Ai(2) is inversely (directly) proportional to /3. Truncated at first K vector 
functional variables a! , . . . , a^- , the p( K ' is said to be the density matrix in the primitive Fourier (PF) approximation. 
As if -4 00, approaches an exact value p at the convergence rate 1/K [341 . 

The structure of p( K ^ is seen to be very similar to that of the classical pdfs and, therefore, we can apply the 
above analytic continuation techniques for thermodynamic averages developed in the classical case. For example, let 
us consider a quantum estimator for the thermodynamic energy which can be obtained from the system partition 
function 



Z {K) (X) = J dxpW(x,x;A). 



(80) 



The expression for the energy is given by 



£/(/?) = - 



dXi 1 8ZW dX 2 1 dZW 



d(3 d/3 ZW dXi d(3 dX 2 

(K + l)d Ai ... 

— ^3 -jHcxoW + McOi(A) 



(81) 
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where the first order cumulants 



Mcio - {S[ K) ) , fJ,coi-\Si K) 

I p,(X,A) \ / p„(X,A) 



exp f-A 1 5f ) (X)-A 2 4 K) (X) 
fdX exp (-XxS A) (X) - A 2 ^ K) (X) 



Here, X labels the whole set of integration variables X =(x, ai, • • • , ajf). Observe that at x = x', s[ K ^ does not depend 
on x. In the case of zero potential energy V = 0, one obtains that /i c xo = cLFT/(2Ai), /x c oi = and U — d/{20). 



Moreover, making use of Theorem lVII.il one can easily derive a quantum estimator for the heat capacity 

~-(K + l)d - 2A 1 // cl0 + \lncw ~ 2*iPn<M + P A*c02 (83) 



k B 

In the limit of zero potential, it is easy to check that \i c \\ = /i C 02 = 0, ^ C 2o = — 9/i c io/9Ai = dK/{2\\) and we 
get Cv/ks = d/2, the value expected for an ideal gas. Similar to expansions (|6"5]l . Theorem IVII.ll can further be 
used to develop Taylor's series expansions of Eqs. (f5Tj) and (|83|l around a fixed value Ao in terms of the higher-order 
cumulants and, thus, to get an analytic continuation of the thermodynamic averages in f3 parameter. 

Note that quantum estimators of the type (|8 1 1) and ((83)) might be more advantageous in MC simulations since 



they require only knowledge of potential functions as opposed to those obtained, e.g., in 27|, |35j, which require first 
and second derivatives of the potential energy to be calculated. It is straightforward with the help of Theorem lVII.il 
to obtain similar analytic continuation formulas in terms of higher order cumulants in the discrctized time-slicing 
primitive approximation [31 1 - 

IX. RESULTS AND DISCUSSION 

To illustrate the above analytic continuation formulas we consider as testing system a cluster of N = 13 neon atoms 
interacting via LJ potential, with the corresponding standard LJ length and energy parameters a^j — 2.749 A and 
clj = 35.6 K being used. The mass of the Ne atom was set to m = 20.0, the rounded atomic mass of the most abundant 
isotope. The particles are assumed to be into the sphere with the confining radius R c = 0.850-LjiV 1 / 3 = 2.0<tlj. 

In Fig. [21 the potential energy pdfs pv(E, T) (labeled by MC) are plotted at fixed temperatures T = 4, 6, 8, 10, 12, 
and 14 K. The MC simulations were imp lemented using the parallel tempering technique, also known as replica 



exchange Markov chain MC sampling j36M44j . In this method several replicas of the same system are simulated in 
parallel in the canonical ensemble, and usually each replica at a different temperature. In this work, 29 replicas have 
been run on the even temperature grid, with the temperature step AT = 1 K, in the interval from 3 to 31 K. Parallel 
tempering is complementary to any set of MC moves for a system at a single temperature, and such single-system 
moves are performed between each attempted swap of complete configurations of the systems at adjacent temperatures. 
The swap moves have been attempted randomly with the probability P SW ap — 0.1/N. The high temperature systems 
are generally able to sample larger volumes of the configuration space, whereas low temperature systems may become 
trapped in local energy minima. Thus, swapping of configurations ensures that the lower temperature systems can 
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access all contributing regions of the configuration integral, thereby overcoming potential barriers between the local 
energy minima. 

The asymptotic formula (O, with the asymptotic parameter r set to 10 2 , has been employed to calculate the 5- 
function. The energy spectra have been calculated on the equidistant grid of 10 3 points in the range [—2000, 1000] 
K. The total number of MC moves has been divided into 50 statistical blocks; the first block data being far from 
equilibration, have been discarded in further averaging. In each block, the number of MC moves has been Nu — 10 5 iV. 
The pdfs in the Gaussian approximation Eq. (|34|) . labeled by Gl, are also shown in Fig. [2j Observe that the Gaussian 
approximation reproduces MC pdfs quite well at higher temperatures T > 12 K, but there both quantitative and 
qualitative differences in the shape of distributions are seen at lower temperatures, especially at T — 8 and 10 K, where 
the "melting peak" in heat capacity Cy is formed. It should be noticed that although 5- function is a non- negative 
function (distribution), its asymptotic approximation (J6j) can in principle be negative at finite values of r. Thus, 
although MC pdfs in Fig [3 are seen to be apparently non- negative, we found that in the regions far from the central 
peaks, where its values are totally corrupted by statistical errors, the MC pdf can take very small negative values 
~ — 10~ 6 — 10~ 8 . In these regions negative MC values should be just zeroed. Also, note that usage of a Gaussian 
representation for the 5- function can be a guarantee for non- negativity of the pdf [suggestion of an anonymous referee] . 

In Fig. [3j the results of inclusion of the 3rd cumulant term Eq. (|37|) . labeled by Airy, are compared with the Gl 
and MC pdfs at temperature To = 8 K. The corresponding MATLAB function has been used to calculate the Airy 
function. In general, the effect of the 3rd cumulant on the pdf results in a slight shift of the peak position to lower 
energies. At To = 8 K one can observe small oscillations in the low energy wing of the Airy pdf. Also, the modulus 
of the Airy curve, labeled by | Airy | , is displayed. 

There is a hint that the MC energy distribution in Fig. [2] at To = 10 K is bimodal. It is believed that such a 
bimodal distribution might have a direct connection to the solid-liquid transition in atomic clusters so that a low- 
energy maximum corresponds to a solid state and a higher one to the liquid state 45 - 4^] - In Fig. |4] (a), the MC 
pdf at To = 10 K is compared to numerical estimates of the integral p v m (|27j) expressed in terms of the cumulant 
expansion (|32[) truncated at k max = 3,5, and 7. One can see that the numerical results including up to the 7th-order 
cumulant term are not sufficient to reproduce the hinted bimodal structure in the MC curve. It is expected that 
inclusion of more cumulant terms will reproduce this structure. Also, numerical results for the statistical temperature 
using the integral representations (|52|) are displayed in Fig. |4] (b). The energy dependence of the temperature is 
seen to be close to a linear one in the neighborhood of fi c i = —1273 K, as qualitatively predicted by Eq. (|56|). Note 
that compution of the statistical temperature at the lower energies E < —1350 K becomes progressively less accurate 
because when moving in the low-energy region far from the central peak, where the pdf values get smaller and, thus, 
relatively less accurate, this ill-defined low-energy region turns out to make a major contribution to the convolution 
integral 

As temperature increases, the difference between the Airy and Gl pdfs becomes negligible. This is demonstrated 
in Fig. [5] at To = 14 K, where Gl and Airy pdfs are compared with the corresponding MC results. 

In Fig. El we demonstrate how the analytic continuation T) — > T formula (|43| works. With the Airy pdf taken 
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at temperature To = 14 K, it is continued to T = 12, 13, 15, and 16 K. For comparison, the corresponding MC 
pdfs are plotted as well. The normalization function has been evaluated with the help of the analytic formula 
(|45p . Obviously, the continuation formula works better when we move to the higher rather than lower temperatures. 
Moreover, the high-energy wings of the continued pdfs are reproduced better than the low-energy ones. Observe that 
the low-energy wing of the Airy pdf at T = 12 K decays faster than the corresponding wing of the MC pdf. The 
increasing difficulty of analytic continuation to the lower temperatures can be explained by the fact that as one can see 
the pdfs are shifted to the lower energies as temperature T goes down so that the overlap between the pdfs at different 
temperatures becomes increasingly smaller. With T decreasing, the difference between inverse temperatures A/3 > 
becomes bigger and the exponential factor in (|43[) grows exponentially. This factor blows up the low-energy wing of 
the p(E, To) pdf and as a result we observe that the peak maximum of p(E, T) gets a shift to lower energies. Moreover, 
in the regions far from the center of the p(E,To) pdf, errors are expected to be dominant. As a result, multiplied by 
an exponentially big factor, these errors either of systematic or statistical nature can generate big deviations from the 
exact values of p(E,T). 

Fig. [7| presents numerical results supporting Theorem IIV.1I First, we calculated moments up to the 7-th order 
on the equidistant grid with the step AT = 1 K in the temperature range from 3 to 31 K. Evaluating higher-order 
moments is computationally inexpensive since it requires only calculation of extra powers of the potential energy. 
Then, with the help of Eqs. (|A8I) cumulants /i c fc, k — 1, . . . , 7 can be recursively obtained. The numerical derivatives 
of cumulants dfj, c k/d/3 are compared with the corresponding higher-order cumulants fi c (k+i) taken with the minus 
sign on the inverse temperature j3 scale. In general, agreement is seen to be better for cumulants of lower orders and 
at higher values of f3. At smaller values of (3 < 0.05 and in the region j3 ~ 0.1 K _1 , where the curves exhibit rapid 
changes, numerical estimates of the derivatives become more scattered due to errors in the numerical formula for the 
derivatives, as well as due to statistical errors present in the MC cumulant estimates themselves. Note that when 
the derivative of fj, c k takes a zero value at some value of [3 ex t, fJ-ck as a function of (3 has an extremum at this point. 
According to Theorem HV.ll fJ. c (k+i) changes the sign at f3 ex t- Fig. [7] confirms such a behavior. 

In Fig. [5] (a)-(c), making use of the Taylor series expansion (|2^|) we present results of analytic continuation of 
the heat capacity calculated at temperatures To = 7, 10 and 14 K, respectively. For comparison, the results of 
MC simulation on the equidistant grid of temperatures with the step AT = 1 K are shown in the range from 3 to 
31 K. The data have been generated with 10 6 iV MC points in each statistical block, and the error bars are at the 
95% confidence level. The present MC data coincide within the statistical errors with the results of independent 
simulations 2jJ obtained in the temperature range 4 to 14 K. Corresponding to the maximum power of A/3 terms 



included in expansion ()26|) . the curves including zero, first and so on up to the 5th-order terms are displayed. One can 
see that dynamics of continuation results is generally improved with inclusion of more terms in the expansion. Thus, 
we find that the results of the 5th-order curve continuation are good in the intervals (6,9), (8.5,11), and (11,31) K 
corresponding to Tq = 7, 10 and 14 K. With increasing temperature To to 14 K, the interval on which the analytic 
continuation works is seen to become bigger. In the neighborhood of To = 10 K, the heat capacity curve achieves 
a maximum value. To find this peak position, we solved the polynomial equations (|25|) for A/3 at /3q — 0.1 K _1 
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truncated at /c max = 1,2,3, and 4. The corresponding roots were obtained using the MATLAB function 'roots'. The 
results are 





= -7.905 
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4 K 




= -7.697 
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= -7.726 
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4 K 




= -7.727 


10" 
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so that the converged peak position temperature is found to be T pea k — 10.078 K. Observe that at To = 7 K the 
slope of the Oth-order curve is negative, while the 1st and higher-order heat capacity curves show positive slopes in 
qualitative agreement with the behavior of the MC curve. 

The Pade approximants often give better approximation of the function than truncating its Taylor series, and it 
may still work where the Taylor series does not converge [48 1 . In Fig. IU(a)-(c), we compare the MC results with the 
Pade approximants of various orders [5/0] , [4/1] , [3/2], and [2/3] calculated in the neighborhoods of temperatures 
To = 7, 10, and 14 K respectively. Notice that the Pade approximant [5/0] coincides with the Taylor series truncated 
at the 5th order. In general, depending on To, the best agreement is seen either for [5/0] or [3, 2] Pade approximants. 
Thus, observe that at To = 10 K the Pade approximant [3/2] shows a slightly better behavior than that of [5/0]. 

In Figs. I10fll2l the corresponding results for the cumulants, their derivatives and the heat capacity curves are 
displayed for N = 55 and 147 neon particles. The confining radius is set to be R c — 0.85<r ljN^^. Notice that the 
MC data have been generated with Nu = 10 6 N MC moves for the system of N = 55 particles and with Nu = 5 • 10 5 A 
MC points for N = 147 particles in each of 50 statistical blocks. The size of the error bars obtained for the cumulants 
in Fig. [10] tends to become bigger for the cumulants of higher orders, but decrease when temperature goes up, except 
for the region where the heat capacity reveals a peak. In Fig. 1 1 1 1 the derivatives d[A c k/ 'd/3, k = 1, . . . , 6 are compared 
to the corresponding cumulants —fi c (k+i)- The agreement is quite good, additionally supporting Theorem II V. 1 1 
Deviations between dfj, c k/df3 and — Mc(fe+i) are m ost pronounced in the regions of cumulant's rapid change and these 
deviations are caused either by the statistical or systematic errors in numerical estimates of the derivatives and the 
corresponding cumulants. Thus, one can see that the 7th-order cumulant for N = 147 particles, shown in the right 
panel of Fig. [T0l is poorly converged at current value of Nu since it is almost totally in error at T < 12 K. As 
expected, this results in a poor agreement observed in the right panel of Fig. II II between the derivative dfi c e/df3 and 
the cumulant —/ic7- Fig. [T^] demonstrates the results of analytic continuation obtained for the heat capacity Cy using 
formula (|26|) . From analytic continuation curves we get the following estimates for the peak positions T pea k = 10.44 
and 12.25 K for the system of N = 55 and 147 particles, respectively. 

Finally, let us consider the convergence issue, that is, how errors present in the MC estimates may affect the results 
of continuation. The MC estimates for cumulants can be represented as /i c /- = ^ e c % act + s c k, where nH act is an exact 
value of the fcth-order cumulant and e c k is an error in its estimate. We assume that cumulant estimates have been 
generated at a fixed number of MC moves. The longer MC moves are generated the less errors one gets in /u c fc's. 
Notice that we cannot directly compare the sizes of errors in cumulants of different orders since they have different 
physical dimensions: [/x c &] = [Energy] fc . In the continuation formula l|26p. the fcth-order cumulant enters multiplied 
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TABLE I: The two standard deviations 2a c k (k — 2, ... ,7) calculated for the fcth-order cumulant at low To = 5, 15 (the size of 
error bars in Fig. I10p and at high temperatures To = 100 and 300 K using 50 statistical blocks. The number of Monte Carlo 
moves in a single statistical block is Nmc — fhicN, where N = 55 is the number of neon particles in the LJ cluster. 
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100 


10 s 


2.5 ■ 


10 2 


1.6 • 10 5 


1.1 • 10 8 


1.2 ■ 


10 11 


1.8 ■ 


10 14 


3.7- 10 17 


100 


10 4 


6.7- 


10 2 


4.7 ■ 10 5 


2.9 ■ 10 8 


2.5 • 


10 11 


2.3 • 


10 14 


2.1 • 10 17 


300 


10 8 


3.9 ■ 


10 2 


6.5 ■ 10 5 


1.7- 10 9 


5.2 • 


10 12 


2.1 • 


10 16 


1.1 • 10 20 



by the fcth power of a small expansion parameter A/3, so that the continued 2nd-order cumulant can be written as 

..exact ( o \ ..exactio \ 

Mc2 QCt (/?) = » e C 2 act (M - Mc3 + ^iM(A/3) 2 + • • • , 

elm = e c2 (M ~ + ^W) 2 + • • • (84) 

One can see that the total error e^W) m the continued 2nd-order cumulant is defined by the sum of errors in the 
2nd- and higher-order cumulants at /3q, contributions from the higher-order cumulant 's errors are being multiplied by 
the powers of the expansion parameter A/3. With the growth of A/3, the contribution of the 3rd- and higher-order 
cumulant's errors to £*2 increases and at some value, which we call a radius of convergence A/3 co „„, its contribution 
becomes comparable to the size of the 2nd-order cumulant's error e C 2{Po)- We do not know exactly errors in cumulants 
(errors are random numbers) but their size can be estimated by 2<7 c fc's, by the two standard deviations in a usual way 
(by error bars in Fig. I10[) . Thus, one can estimate the raduis of convergence Aftcorn?^ for the fcth-order cumulant 
(fc > 2) by the expression 

/ rr \ V(*- 2 ) 

|A/3(^)|=((fe-2)!^ (85) 

Using the relationship between inverse and direct temperature scales A/3 = j3 - /3 = 1/T - 1/T « -AT/T 2 , where 
AT = T — T , one finds the corresponding radius of convergence in the temperature scale 

, x l/(fc-2) 

ATi^^T 2 (fc-2)!^ (86) 

V a ck) 

It roughly defines the temperature range within which the errors present in the fcth-order cumulant will produce the 
same order errors as in the ii c iifio)- 

The 2cr c fc's obtained for the cluster of N — 55 neon particles at low To = 5 and 15 K and high temperatures 
T = 100 and 300 K are summarized in Table HI Observe how the magnitude of errors in cumulants increases, roughly 
by an order, as the number of MC moves specified by the parameter Jmc decreases by two orders at T = 100 K. 

Making use of the data in Table HI we can estimate radii of convergence for the higher-order (3-7) cumulants with 
the help of Eq. (|56"|) ; the results of evaluation are summarized in Table [TT1 Observe that at T = 5 K, the higher-order, 
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TABLE II: The radius of convergence for the fcth-order cumulant Ario^J 2 ' calculated with the help of Eq. (|86[) an d the data 
from Table [1] Other notations are the same as in Tabled] 



To (K) 


!mc 


AT^ V 2) (K) 


ATi 4 ^ 2) (K) 


ATi 5 ^ 2 ' (K) 


ATi 6 ^ 2) (K) 


L-*-L conv 


5 


10 6 


0.3 


0.4 


0.6 


0.06 


0.04 


15 


10 6 


0.8 


0.9 


1.1 


1.2 


1.2 


100 


10 6 


25 


24 


26 


29 


30 


100 


10 5 


16 


21 


23 


24 


24 


100 


10 4 


14 


21 


25 


29 


33 


300 


10 6 


54 


61 


69 


74 


76 



6th and 7th-order, cumulants appear to be non-converged since their radii of convergence are by an order of magnitude 
smaller than the values obtained for the 3rd-, 4th-, and 5th-order cumulants. At higher temperatures shown in Table 
III1 all the cumulants seem to be well converged; the smallest radius of convergence is generally found for the 3rd-order 
cumulant 's contribution. As an empirical rule, we find that the radius of convergence is roughly proportional to the 
temperature To. Notice that at Tq = 100 K, when the number of MC moves drops down by two orders, the radius of 
convergence ATion^ 2 " 1 is reduced by about two times only. 



X. CONCLUDING REMARKS 



In obtaining reliable estimates of thermal averages in realistic, multidimensional, classical or quantum many- or 
few-body systems, the MC simulation is often the only way of getting right answers. The simulated averages depend on 
thermodynamic parameters, such as temperature, volume etc., and/or interaction parameters between the particles. 
Therefore, in order to avoid many time-consuming runs of the MC codes at various parameters, it is important to 
develop robust analytic techniques that allow us to continue the MC data in the neighborhood of a set of prescribed 
parameter values. The key finding of this work is a simple relationship Eq. (|21|) between derivatives and the higher- 
order cumulants. Since many important thermodynamic quantities, such as energy and heat capacity, can be expressed 
in terms of cumulants, this theorem provides an analytic tool or bridge to continue MC data in the neighborhood of 
a control parameter value, say, /3rj or to fill a gap by constructing an analytic bridge in between of two neighboring 
parameters /3q and j3±. 

To find an optimal numerical scheme to evaluate the cumulants up to the fc maa; th-order (in our examples k max = 7), 
it is important to estimate the amount of numerical work required. A reasonable estimate of this time is the number 
of potential function calls N ca u required to compute a thermodynamic quantity, e.g., the heat capacity. In a single 
MC move, the system moves to a new position and one has to calculate the potential function Vi = V(Xi) one time 
at a new position X, so that N ca u = Nu where Nu is a total number of MC moves in a single statistical block. To 
calculate cumulants up to the fc moa: th-order, one has to compute additionally k max — 1 powers Vj° , k = 2, . . . , k max or 
perform k max — 1 multiplications of the known potential function value Vi . These multiplications are computationally 
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cheap and do not depend on the size of the system. Therefore, as estimated the amount of numerical work to evaluate 
the higher-order cumulants will be practically the same as in the standard parallel tempering scheme which requires 
evaluation of only Vi and values. 

By Eqs. (j2"T)) - ([3"Tj) . the energy pdfs can also be expressed in terms of cumulants. Truncating the cumulant expansion 
at k max , we have derived analytic expressions ()33j) . (|34j) . and (|37|) for the pdfs at k max = 1,2, and 3 respectively. The 
higher-order energy pdfs truncated at k max > 4 can be obtained either by evaluating the r-integral numerically or by 
applying an appropriate (saddle-point) asymptotic method. 

Generalizations to a multi-parameter classic system are rather straightforward and can be expressed by the Theorem 
IVII.ll Using this theorem, one can easily develop similar analytic continuation formulas, such as Eqs. (|68[) . in 
terms of the higher-order multi-variate cumulants in the multi-parametric space. The path integral Feynman-Kac 
representation of the density matrix is a very useful formulation of the quantum statistical equilibrium operator since 
it basically reduces the problem of analytic continuation in quantum case to the corresponding problem in multi- 
parametric classical systems. Of course, technically this reduction can be done if the original infinite-dimensional 
integration over the path variables can be replaced somehow by a finite-dimensional one. We considered several 
possible scenarios, having different convergence properties with respect to the number of path variables, of such 
infinite-to-finite dimensional integral replacements in Section [Villi 

Numerical testing of Theorem I IV. 11 analytic continuation formulas for the energy pdfs and heat capacity curves 
has been exemplified by an LJ classical cluster sytem in Section IIXI By now, usefulness of analytic continuation 
formulas in multi-parametric classical and quantum systems have not yet been tested numerically. Further numerical 
investigations of interesting multi-parameter systems are of special interest. The results of such investigations when 
available will be published elsewhere. 



This work was supported by NRF (National Honor Scientist Program: 2010-0020414, WCU: R32-2008-000-10180-0) 
and KISTI (KSC-2011-G3-02). 



Relationships between univariate cumulants and moments can be derived from the moment-generating function [1 
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Appendix A: Moments versus Cumulants 



M(t) = (exp(tfl)) = 1 + 




(Al) 



where [if. = (i? fc ) are called moments. Assuming the average of 1 to be nonzero, we conclude that In (exp (tH)) is an 
analytic function of t, in a vicinity of zero. Therefore, it has a Taylor expansion with respect to t. This expansion is 
called cumulant expansion; the coefficients [i c k of the expansion are called cumulants. Thus, cumulants in terms of 
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moments are expressible by the formula 

fc 



(-l) r+1 ^ fc! 

Mcfe = 2^ — - — 2^ „.■...„ . /v •••/'/• ( A2 ) 

r=l ' pi,...,p r GC fer 



where the summation of the integer indices pi , • • • } p r > 1 is restricted by the relation 

Cfer = j(pi, ' ' ' ,Pr) ■ Y^Pj = k 

In order to calculate the fc-th order cumulant one needs moments up to the fc-th order. 
On the other hand, moments can be expressed via cumulants by a similar formula 

1 fc! 

«=Eh E P1 i..'. p i fa-^ ( A3 ) 

r=l ' pi,..., Pr ec tr Fi ' 

The fcth moment is a fcth-degree polynomial in the first fc cumulants. These polynomials have a remarkable 
combinatorial interpretation: the coefficients count certain partitions of sets. A general form of these polynomials is 

^ k = X II Ms\ (A4) 

7T SG7T 

where 7r runs through the list of all partitions of a set of size fc; "S 1 € 7r" means 5 is one of the "blocks" into which 
the set is partitioned; and \S\ is the size of the set 5. To better understand how Eq. (|A4j) corresponds to (|A3| let us 
consider all possible partitions of the set of natural numbers {1, 2, . . . , fc}. We can divide all possible partitions {tt} 
into disjoint groups {n} = {tti, . . . , ~Kk\ such that the number of blocks in a particular partition 7iy, r = 1, . . . , k is 
equal to r. For example, in case of fc = 3 we obtain the set of all possible partitions classified as {tti, 7T2, 7^}. Here, 7Ti 
includes an "improper" partition (123). The size of the (123) partition is three. 7T2 corresponds to partitions into two 
blocks: (12)3, (13)2, (23)1, with sizes of blocks being two and one. tt^ splits the set into three blocks (1)(2)(3), each 
block of size one. Round parenthesis show how the set is partitioned. The total number of partitions of a fc-element 
set is the Bell number Bk'- Bq = 1, B\ = 1, Bi = 2, B3 = 5. Bell numbers satisfy the recursion Bk+i = 53 ra =o (fi)^"' 
where ( ) is the binomial coefficient. 

If Si , . . . , S r € 7r r are r blocks into which the set can be divided, then Eq. (|A4[) can be rewritten as 

fc 

A*fc = X X MSi\---MSr\ ( A5 ) 

r=l Si,...,S r er r 

where inner summation runs over all possible non-empty, disjoint blocks. Then, each r-term in Eq. (IA3|) corresponds 
to the r-term in (|A5|) : 

fc! 



^ X „J,,,„J ^P1 •••P-CPr - X ^C|S!| •••Mc|S r | (A6) 

' Pi, •••,?>. 

Indeed, the multinomial coefficient 



r ' pi,...,p r £C kr Pl ' "' Pr ' Si,...,s r e7iy 
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in the left-hand-side is the number of ways of grouping k objects (numbers) into r groups (blocks) of sizes pi, . . . ,p r , 
when the order within each group does not matter. The summation over blocks in the right-hand side can be carried 
out into two steps. First, the summation over all possible block's sizes \Si\ = p%, . . . , \S r \ — p r such that p\ < • • • < p r 
and then summing up over blocks at fixed sizes of blocks 

E •••= E E - 

S x ,...,S T £lt r Pl<---<Pr |Si|=pi,...,|S r |=p r 
plH hpr = fc 

Obviously, the inner summation will result in the same multinomial coefficient (|A7j) as in the left-hand side of (|A6|) , 
Finally, notice that the summed function fi cpi ■ ■ ■ fi cPr , and the multinomial coefficient are totally symmetric functions 
with respect to permutations of p\, . . . ,p r indexes and r! is the number of their permutations. This symmetry results 
in 

E ■ • ■ — E 

Pl>- ,Pr Pl<-<Pr 
PlH hPr=fc pi _| |-Pr=fc 

This proves Eq. (|A"6|) . 

Using Eq. (|A3I) or (IA5|) . one obtains the following expressions for the first seven moments in terms of cumulants 
A'i = Mcl) 

M3 = Mc3 + i^c2[lc\ + Mel! 

= Mc4 + 4/i c3 Ai c l + 3^2 + 6A*c2A*?l + Mel) 
^5 = /U c5 + 5/i c 4Aicl + 10^ c 2Mc3 + 10^c3Aicl + 15 M?2Mcl + 10/ic2Mcl + Mcli 
M6 = Mc6 + 6/i c5 /i c i + 15^ c2 ^c4 + 10^c 3 + 15Aic4Mcl 

+ 60/i c3 /i c2 Ai c i + 15^2 + 2 Q/"c3A*ci + 45/4 2 Mci + McU 
M7 = Mc7 + 7/i c6 ^i c i + 21^ c2 ^c5 + 35^ c3 ^c4 + 2f/x c5/ u^ 1 

+ f 05/! c 4Ai c 2Mel + 70/ic3Mcl + 105^c3/l^ 2 + 35^4^^! 

+ 210/i c2 Ai c3 Mci + 105/ic2Mci + 35m c3 Mc1 + 105Mc2Mci + 21/^/4 + Mcl ( A $) 

These equations relating the cumulants and the moments can also be interpreted as recurrence relations that allow 
the expression of the higher-order cumulants in terms of lower-order ones. 

Appendix B: Proof of Eq. (pTTj) 

The proof is by induction. At k = 1, the statement of theorem follows directly from Eq. (|2TJ)) and definitions of 
the first two cumulants /i c i = /ii and \i C 2 = /i2 — Mi- Let us assume that Eq. (|2ip is valid at k = 1, . . . ,p; we have 
to prove its validity at k = p + 1. Making use of combinatorial representation (|A5[) . the (p + I )-st moment can be 
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wriiicn as 1+ 

P+i 

Mp+i = E E MciSii •••Mc|S r | (Bl) 

r=l S lt ...,S T 

where S\, . . . ,S r denote a partition of a set of natural numbers {1, . . . ,p + 1}. If blocks S\, . . . ,S r contain, respectively, 
fci, . . . , k r elements [numbers] such that k\ + . . . + k r = p + 1, then /J, c \Si\ — A*cfei > ■ • ■ ; Mc|s r = ^ck T - The term in (|Blj) 
for r = 1 corresponds to /i c ( P +i), so that Eq. (|B1[) can be rewritten as 

p+i 

Mc(p+1) = Mp+1 ~ E E •••Mc|5r| ( B2 ) 

r=2 Si,...,S r 

where the summation is performed over the "proper partitions" when r > 2. Observe that the sizes ki,...,k r of 
"proper partitions" in (|B2|) cannot be bigger than p. Differentiating (|B2[) with respect to /3, we get 

^Mcfp+l) 



p+1 



Mp+2 + Mp+l/^l 



r=2Si,...,5 r K 

-E E /^i-^f (B3) 

where in the first line we have used Eq. (1201) at k = p + 1. The term /Lt p +i/ui can be rewritten as 

p+i 

/^ + l/!l=^ ^ Mc|S 1 |---Mc|Sr.|Mc|T| (B4) 
r=l Si,...,S r 

Here, to the set of natural numbers {1, . . . ,p + 1} we added the number p + 2 so that in the partition Si, . . . , S r , T 
the 5"s correspond to all possible blocks partitioning the set {1, . . . ,p + 1} and T is related to a fixed element, the 
number p + 2. Evidently, \T\ = 1 and /x c m = /x c i = /ii- Further, by induction conjecture one can use Eq. (|2~T|) for 
the cumulant derivatives in the second and next lines 

MciSii 

— ~ Mc(fci+1) — -A*c|SiUT|! 

Me [Sr] 



-Mc(fc r +i) — -^c\s r ur\, (B5) 



d/3 

where U T, . . . , S r U T denote that to a particular partition of {1, . . . ,p + 1} we add consecutively a fixed element 
T = p + 2 to Si, . ■ . , S r blocks. Substituting Eqs. $B4$ and (|B5|) into fB3]). we get 



, p+i 
«A*c(p+l) , 

To = -Mp+2 + 2^ MSi\ • ■ • t*c\Sr\Vc\T\ 

P r=lSi,...,S r 

P+1 

+ E E Mc|SlUT| ■ --MclSrl H 

r=2 Si,...,S r 
p+1 

+ E E MclSil •••Mc|S r uT| (B6) 

r=2 Si,...,S r 
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Evidently, summing over Si, . . . , S r partitions of{l,...,p+l} with inclusion of a fixed element T is equivalent to 

p+2 

E E 

r=2 S lt ...,S r 

where Si, . . . ,S r is a partition of {1, . . . ,p + 2} at r > 2, so that Eq. (|B6j) can be rewritten as 

p+2 



d/3 



"Mp+2 + E E Mc|Si| •••Mc|S,| 
r=2 Si,...,S r 

-Mc(p+2) 



(B7) 



Here, in transition to the second line, we employed Eq. (|B2j) . where a substitution p — > p + 1 has been made. This 
completes the proof. 



Appendix C: Multivariate Cumulants and Proof of Theorem I VII. H 



Similar to Eq. (|Al|) . multivariate moments and cumulants are defined via the moment generating function 

/ / P \ \ ,ki t k P 

M(ti,...,t p ) = (exp [^t k H k ) = !+ J! <■ .-■ 



\k=l 



/.•,!•••/.•,.!' 



exp 



ifci <fep 

Et-j • • • hp 



fcl ,. . . , 



where fJ,ki,...,k p — \ H± 1 • • • ) and summation runs over non-negative integers ki, . • . , k p except for k\ 



(Cl) 



0. Using Eq. (|C1|) . one can derive explicit relations between p-dimensional moments and cumulants, multivariate 
generalizations of Eqs. (IA1|) and (|A3I) ; see Appendix in [29(. Following 29(, it is convenient to introduce compact 
notations for multi-indexes k = {ki, ■ ■ . , k p ), as well as for the product of powers of multivariate variables t k = 
ij 1 • • • tp F , and product of factorials k! = kil ■ ■ ■ k p \. In these notations, Eq. (|CJ1[) reads 



y- 

4^ k! 



/i k = exp 



En' 1 * 



E 1 E 1 fMciM ' • • Mc 

r i m! • • • n r ! 



(C2) 



r=l ' ni^O,...,n r ^O 

In the second line, the exponential function has been expanded in Taylor's series and multi-indexes in the inner-sum 
are non-negative such that ni ^ 0, ...,n r ^ 0. Differentiation of Eq. (|C2[) by applying the multi-index differential 
operator 

g|k| 

(C3) 



Dt 

dtl 1 ■ ■ ■ dtp" 

to the both sides and then putting t = yields a multivariate analogue of (|A3 

|k| 



= E ^ E 

r=l m^O,...,n r ^O 



k! 

ni! • • • n r ! 



P'Clli ' ' ' P'C 



(C4) 
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where |k| = k\ + . . . + k p . 

In multi-dimensional case, Eq. (|A5j) remains valid if the single index k is replaced by a multi-dimensional 
one k, and partitioning 7iy of the set {1,2, ... ,k} is replaced by partitioning of the ordered set of multi- indexes 
({1, . . . , fci}; {1, - . . , fc 2 }; ...;{!,..., fc p }): 

|k| 

Mk = ^2 E ' "MSr\ 



r = l S 1 ,...,S r eiTr 

|k| 

EE E 

r=l m<---<n r |Si|=m,...,|S r |=n r 
niH hn r =k 



(C5) 



where S\,...,S r are all possible blocks into which the set of multi- indexes can be divided. The summation over 
blocks in the second line is organized as follows. The first summation over r defines the number of blocks S\, . . . ,S r 
partitioning the set of multi-indexes. The second summation runs over all possible sizes of blocks such that ni + 
• • • + n r = k. In components, if = (rin, Uj2, • • ■ , nip), i — 1, . . . , r [the first subindex labels the block's number, the 
second one is due to variable's number], we have constraining equations 



ni 



n r = k 



nn + • • • + n r i = ki 
?ii2 + • • • + n r2 = k 2 



(C6) 



and inequalities specifying the partial ordering conditions 



ni < • • • < n r 



ni2 < n 2 2 < ■ ■ ■ < n r2 



(C7) 



ni p < n 2p < 

The inner summation over all possible blocks at fixed sizes results in 

k! 



|Si|=ni,...,|S r .|=n r 



ni 



!...n„l 



< n r 



(C8) 



where 



k! 



k x \---k p \ 



ni! • • • n r ! n\\\ ■ ■ ■ n\ p \ ■ ■ ■ n r \\ ■ ■ ■ ,L Tp . 

Observe that the multi- index- multinomial coefficient (|C9j) in the right-hand side of (|C8|) is the same as in (|C4j) . As in 
Appendix XF\ due to symmetry of the summed functions with respect to permutations of multi-indexes ni, . . . , n r , we 
obtain that Eqs. (|C4[) and (|C5[) are equivalent representations for multivariate moments. The latter representation is 
helpful in proving the multivariate analogue of the theorem on cumulant's derivative. 



,!..., 



I . . . . 



(C9) 
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The proof of Theorem IVII.ll is by induction and quite similar to the proof of Theorem IIV.1I done in previous 
Appendix. The only complication is due to multi-index nomenclature. 
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FIG. 2: The potential energy pdfs pv(E, T) for the system of TV = 13 neon particles with atomic mass m — 20.0 a.u. interacting 
via Lennard- Jones potential with the parameters ol.i = 2.749 A and £lj = 35.6 K as a function of energy E calculated at 
fixed temperatures T = 4, 6, 8, 10, 12, and 14 K. The particles are assumed to be in the sphere with the confining radius 
R c = 2.0(j_lj. MC labels the total results of MC simulations using the asymptotic formula ([6]) at r = 10 2 for S- function, whereas 
Gl is the Gaussian approximation Eq. (|34[) to the MC pdf. 
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FIG. 3: Effect of the 3rd cumulant term on the pdf at T = 8 K. Here, 'Airy' labels the Airy pdf Eq. (l37l) . whereas ' | Airy | ' is 
the modulus of the Airy pdf. The system and other notations are the same as in Fig. [2] 
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FIG. 4: Effects of the higher-order, k max = 3, 5, and 7 cumulants on the pdf (a) and on the statistical temperature (b) at 
To = 10 K. Here, pj, and T^J^y*^ are the pdfs and the statistical temperatures, respectively, calculated with up to the 
fcmaK-order cumulants included. The MC pdf shows some bimodal structure. The system and other notations are the same as 
in Fig. H 
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FIG. 5: The Airy and the Gaussian pdfs in comparison with the MC pdf at To = 14 K. The third cumulant term has a minor 
effect on the pdf at To = 14 K. The system and notations are the same as in Fig. [3] 
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FIG. 6: The analytic results from formula (I43p applied to the Airy pdf at To = 14 and continued to temperatures T = 12, 13, 15, 
and 16 K in comparison with the corresponding MC results. The system and notations are the same as in Fig. [3] 
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, 6. Numerical estimates of the derivatives dfick /d/3 
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FIG. 8: The heat capacity Cy in units of the Boltzmann constant fcs as a function of temperature T. The system is the same 
as in Fig. [2] The data labeled by MC have been generated with l(f N MC points in each of 50 statistical blocks, and the error 
bars are at 95% confidence level. The analytic continuation of the heat capacity value in the neighborhoods of To = 7 (a), 10 
(b), and 14 K (c) are done using expansion formula (|26[) . The fcth-order curves (k = 0, . . . , 5) are the results of calculation that 
include maximum up to the kth power terms in the expansion ()26[) in A/3 powers. 
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FIG. 9: Analytic continuation of the heat capacity by Pade approximants of the orders [5/0], [4,1], [3,2], and [2,3] in the 
neighborhoods of To = 7 (a), 10 (b), and 14 K (c). The system and other labels are the same as in Fig. [8] 
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FIG. 10: The temperature dependence of cumulants /x c fc, k = 1, . . . , 7, calculated for the system of iV = 55 (left) and 147 (right 
panel) neon particles. The MC data have been generated with 10 6 iV and 5 ■ 10 5 iV MC points in the left and right panels, 
respectively, in each of 50 statistical blocks. The error bars are at 95% confidence level. 
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FIG. 11: Numerical estimates of the derivatives d^, c k/d/3, k = 1, . . . , 6 (red circles) as compared to — ^i c ^+i) (black crosses). 
The system and statistics are the same as in Fig. 1101 
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FIG. 12: The results of analytic continuation obtained with the help of formula (|26[) for the system of iV = 55 (upper) and 147 
(lower panel) neon particles. Notations are the same as in Fig. [8] Statistics of the MC moves are the same as in Fig. 1101 



